logo_datactivist           logo_pleiade


Objectif : Exploration des données et modélisations

Date début de l’analyse : 17 mai 2022

Date de la dernière modification : 14 juin 2022

# Packages
packages = c("tidyverse", "glue", "kableExtra", "knitr", "plotly", "hrbrthemes", "gridExtra", "stats", "arrow", "lme4", "tidymodels", "broom.mixed", "multilevelmod", "dotwhisker", "rAmCharts", "DescTools")

## Installation des packages si besoin et chargement des librairies
package.check <- lapply(packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)}})

# Import BSO complet (2013-2020)
data <- read_csv("../data/4.process/managing-variables_2nd-part/whole_bso.csv", col_types = cols(author_position = "c"), na = c("", "NA"))

I - Modélisations for french AND non french CA articles


Analyse exploratoire


Péparation des données aux modélisations

# Création de nouvelles variables : couleur de journal finale et APC fiables
data <- data %>% mutate(oa_color_finale = case_when(oa_color_article_BSO == "diamond" ~ "diamond",
                                                    oa_color_openalex == "gold" ~ "gold",
                                                    oa_color_openalex == "hybrid" ~ "hybridOA",
                                                    TRUE ~ "other"), 
                        amount_apc_EUR_trusted = case_when(apc_source == "openAPC_estimation_publisher" | apc_source == "openAPC_estimation_publisher_year" ~ NA_real_,
                                                           TRUE ~ amount_apc_EUR))

Pour les modélisations, nous cherchons à expliquer et prédire le montant des APC (Y) avec 2 variables potentielles :

  • amount_apc_EUR : variable du BSO sans modification ;
  • amount_apc_EUR_trusted : montant des APC fiable, c’est-à-dire dont la source est doaj, OpenAPC, openAPC_estimation_issn_year ou openAPC_estimation_issn. Sur les 6 sources d’APC du BSO, “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year” sont écartées.

L’analyse exploratoire nous dira quelle variable des APC est à retenir pour inclure dans les modèles. Dans tous les cas, les variables explicatives utilisées pour prédire Y seront “oa_color_finale”, “tier”, “bso_classification”, “is_french_CA_bso_wos_openalex_single_lang”, “year” et “is_covered_by_couperin”.

Pour préparer les données nous créons donc les variables “oa_color_finale” et “amount_apc_EUR_trusted”, puis nous sélectionnons les 6 variables qui entrent dans les modèles en supprimant les valeurs manquantes, et enfin nous filtrons sur les articles pour lesquels des APC ont été payés (apc_has_been_paid == 1, amount_apc_EUR_trusted > 0).

# Sélection de variables pour les modèles
bso_mvp <- data %>% 
    select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR, amount_apc_EUR_trusted, is_covered_by_couperin) %>% 
    mutate(bso_classification = as.factor(bso_classification), 
           apc_has_been_paid = as.factor(apc_has_been_paid), 
           tier = as.factor(tier), 
           oa_color_finale = as.factor(oa_color_finale), 
           # On modifie les valeurs des catégories pour éviter les conflits de nom de classes dans les modélisations
           is_french_CA_bso_wos_openalex_single_lang = as.factor(case_when(is_french_CA_bso_wos_openalex_single_lang == 1 ~ "french CA",
                                                                           is_french_CA_bso_wos_openalex_single_lang == 0 ~ "non french CA")), 
           is_covered_by_couperin = as.factor(case_when(is_covered_by_couperin == 1 ~ "covered",
                                                        is_covered_by_couperin == 0 ~ "non covered")))

# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>% 
    filter(!is.na(bso_classification),
           !is.na(oa_color_finale),
           !is.na(is_french_CA_bso_wos_openalex_single_lang),
           !is.na(amount_apc_EUR),
           !is.na(is_covered_by_couperin)) %>% 
    filter(apc_has_been_paid == 1, amount_apc_EUR > 0)


Y moyen par année et par variable catégorielle


Y1 = amount_apc_EUR

Discipline
table_y1_class <- bso_pop %>% group_by(year, bso_classification) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR)) %>% ungroup()

  # Plot
graph_y1_class <- ggplot(table_y1_class, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y1_class)


Tier
table_y1_tier <- bso_pop %>% group_by(year, tier) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR)) %>% ungroup() %>% mutate(tier = as.character(tier))

  # Plot
graph_y1_tier <- ggplot(table_y1_tier, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y1_tier)


French CA
table_y1_ca <- bso_pop %>% 
    group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>% 
    summarise(montant_moyen = mean(amount_apc_EUR)) %>% ungroup() %>% 
    mutate(is_french_CA = as.character(is_french_CA_bso_wos_openalex_single_lang))

  # Plot
graph_y1_ca <- ggplot(table_y1_ca, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y1_ca)


Couleur OA
table_y1_color <- bso_pop %>% group_by(year, oa_color_finale) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR)) %>% ungroup() %>% mutate(oa_color_finale = as.character(oa_color_finale))

  # Plot
graph_y1_color <- ggplot(table_y1_color, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y1_color)


Covered by Couperin
table_y1_coup <- bso_pop %>% group_by(year, is_covered_by_couperin) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))

  # Plot
graph_y1_coup <- ggplot(table_y1_coup, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y1_coup)


Y2 = amount_apc_EUR_trusted

Discipline
table_y2_class <- bso_pop %>%  group_by(year, bso_classification) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()

  # Plot
graph_y2_class <- ggplot(table_y2_class, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y2_class)


Tier
table_y2_tier <- bso_pop %>% group_by(year, tier) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(tier = as.character(tier))

  # Plot
graph_y2_tier <- ggplot(table_y2_tier, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y2_tier)


French CA
table_y2_ca <- bso_pop %>% 
    group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>% 
    summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% 
    mutate(is_french_CA = as.character(is_french_CA_bso_wos_openalex_single_lang))

  # Plot
graph_y2_ca <- ggplot(table_y2_ca, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y2_ca)


Couleur OA
table_y2_color <- bso_pop %>% group_by(year, oa_color_finale) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(oa_color_finale = as.character(oa_color_finale))

  # Plot
graph_y2_color <- ggplot(table_y2_color, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y2_color)


Covered by Couperin
table_y2_coup <- bso_pop %>% group_by(year, is_covered_by_couperin) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))

  # Plot
graph_y2_coup <- ggplot(table_y2_coup, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y2_coup)


Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable. Nous supprimons donc la variable des APC du BSO puis filtrons sur les articles ayant un APC fiable connu.

n_old <- nrow(bso_pop)
bso_pop <- bso_pop %>% select(-amount_apc_EUR) %>% filter(!is.na(amount_apc_EUR_trusted))

Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et des valeurs connues pour les autres variables.


Valeurs extrèmes

amBoxplot(bso_pop$amount_apc_EUR_trusted, xlab=" ", ylab="Montant des APC fiables", main="Boxplot de Y")

L’étendue (écart entre les valeurs minimales et maximales) du montant des APC est importante, avec des valeurs allant de 2.27 à 9848.4 euros. Nous constatons 2 valeurs extrêmes qui se détachent largement des autres sur le boxplot, où les montants sont respectivement de 9848.4 et 9028.43 euros. Après vérification à partir de l’ISSN, les valeurs sont plausibles donc nous décidons de laisser ces valeurs aux données à modéliser.


Normalité de distribution

# Normalité de Y
mean <- mean(bso_pop$amount_apc_EUR_trusted, na.rm = TRUE)   
sd <- sd(bso_pop$amount_apc_EUR_trusted, na.rm = TRUE)  
ggplot(bso_pop, aes(bso_pop$amount_apc_EUR_trusted)) +
      geom_histogram(aes(y=..density..),  color="black", fill = "steelblue", alpha = 0.2) +
    geom_density( color="red", size = 1, adjust = 5) +
    stat_function(fun = dnorm, colour = "Black", size = 1, args = list(mean = mean, sd = sd))  +
    xlim(c(min(bso_pop$amount_apc_EUR_trusted)-1, max(bso_pop$amount_apc_EUR_trusted)+1))+ ggtitle("Distribution du montant des APC fiables") +
  xlab("Montant des APC fiables") + ylab("Densité") +
     theme_classic()

# Test statistique de normalité
out <- JarqueBeraTest(bso_pop$amount_apc_EUR_trusted, robust = FALSE, method = "chisq")  # Y non normal (on rejette H0)
library(formattable)
df <- data.frame(value = unlist(out))
tdf <- as.data.frame(t(df))
formattable(tdf)
statistic.X-squared parameter.df p.value method data.name
value 49797.6693637252 2 0 Jarque Bera Test bso_pop$amount_apc_EUR_trusted

Le montant des APC fiables n’est pas normalement distribué.


Corrélations et dépendances entre les variables

# Fonction pour sauvegarder les résultats des tests
save_results <- function(output, num){
    df <- data.frame(value = unlist(output))
    tdf <- as.data.frame(t(df))
    assign(glue("tdf_{num}"), tdf[,1:5], envir = .GlobalEnv)
}
# Test de corrélation (Spearman car distribution non normale)
output <- cor.test(bso_pop$amount_apc_EUR_trusted, bso_pop$year, method = "spearman") # corrélation significative
df <- data.frame(value = unlist(output))
tdf <- as.data.frame(t(df))
tdf <- tdf %>% mutate(data.name = str_replace_all(data.name, "_", "\\_"),  #escape special characters
                          data.name = gsub("bso_pop\\$", "", data.name))
formattable(tdf)
statistic.S p.value estimate.rho null.value.rho alternative method data.name
value 242633003532012 0 0.149384487346854 0 two.sided Spearman’s rank correlation rho amount_apc_EUR_trusted and year 
# Tests de dépendance
    # Année
save_results(chisq.test(bso_pop$year, bso_pop$bso_classification), "1")
save_results(chisq.test(bso_pop$year, bso_pop$tier), "2")
save_results(chisq.test(bso_pop$year, bso_pop$oa_color_finale), "3")
save_results(chisq.test(bso_pop$year, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "4")
save_results(chisq.test(bso_pop$year, bso_pop$is_covered_by_couperin), "5")
    # Discipline
save_results(chisq.test(bso_pop$bso_classification, bso_pop$tier), "6")
save_results(chisq.test(bso_pop$bso_classification, bso_pop$oa_color_finale), "7")
save_results(chisq.test(bso_pop$bso_classification, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "8")
save_results(chisq.test(bso_pop$bso_classification, bso_pop$is_covered_by_couperin), "9")
    # Tier
save_results(chisq.test(bso_pop$tier, bso_pop$oa_color_finale), "10")
save_results(chisq.test(bso_pop$tier, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "11")
save_results(chisq.test(bso_pop$tier, bso_pop$is_covered_by_couperin), "12")
    # Couleur OA
save_results(chisq.test(bso_pop$oa_color_finale, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "13")
save_results(chisq.test(bso_pop$oa_color_finale, bso_pop$is_covered_by_couperin), "14")
    # French CA
save_results(chisq.test(bso_pop$is_french_CA_bso_wos_openalex_single_lang, bso_pop$is_covered_by_couperin), "15") #toutes dépendantes
# Présentation des résultats dans une table
output <- rbind(tdf_1,tdf_2,tdf_3,tdf_4,tdf_5,tdf_6,tdf_7,tdf_8,tdf_9,tdf_10,tdf_11,tdf_12,tdf_13,tdf_14,tdf_15) %>% 
    mutate(p.value = as.numeric(p.value),
           Dependance = case_when(p.value <= 0.05 ~ "Oui",
                                  TRUE ~ "Non"))
# Escape special characters
output <- output %>% mutate(data.name = str_replace_all(data.name, "_", "\\_"),  #escape special characters
                          data.name = gsub("bso_pop\\$", "", data.name))  #remove base name
rownames(output) <- NULL  #remove rownames
# Display table
kable(output, format = "html", caption = "Résultats des tests de dépendance entre les variables qualitatives") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
Résultats des tests de dépendance entre les variables qualitatives
statistic.X-squared parameter.df p.value method data.name Dependance
1844.04364295786 63 0 Pearson’s Chi-squared test year and bso_classification Oui
10564.3561750721 21 0 Pearson’s Chi-squared test year and tier Oui
89.5375126042915 7 0 Pearson’s Chi-squared test year and oa_color_finale Oui
207.441377682062 7 0 Pearson’s Chi-squared test year and is_french_CA_bso_wos_openalex_single_lang Oui
155.842642000071 7 0 Pearson’s Chi-squared test year and is_covered_by_couperin Oui
7903.46367919093 27 0 Pearson’s Chi-squared test bso_classification and tier Oui
4498.17074461241 9 0 Pearson’s Chi-squared test bso_classification and oa_color_finale Oui
223.795592556606 9 0 Pearson’s Chi-squared test bso_classification and is_french_CA_bso_wos_openalex_single_lang Oui
13109.2394339774 9 0 Pearson’s Chi-squared test bso_classification and is_covered_by_couperin Oui
11883.1772351368 3 0 Pearson’s Chi-squared test tier and oa_color_finale Oui
228.437752033453 3 0 Pearson’s Chi-squared test tier and is_french_CA_bso_wos_openalex_single_lang Oui
15863.9058727782 3 0 Pearson’s Chi-squared test tier and is_covered_by_couperin Oui
1870.71182474937 1 0 Pearson’s Chi-squared test with Yates’ continuity correction oa_color_finale and is_french_CA_bso_wos_openalex_single_lang Oui
49839.9633692406 1 0 Pearson’s Chi-squared test with Yates’ continuity correction oa_color_finale and is_covered_by_couperin Oui
1145.55413773546 1 0 Pearson’s Chi-squared test with Yates’ continuity correction is_french_CA_bso_wos_openalex_single_lang and is_covered_by_couperin Oui 
# Vérifier : 
    #- homoscédasticité (homogénéité de la variance au sein des sous-pops)
    #- indépendance (sous-pops doivent avoir un effet aléatoire sur Y)
    #- autocorrélation (résidus et erreurs ne doivent pas être auto-correlés)


Modélisations


Régressions linéaires : effets fixes


Modèle 1

  • Y : amount_apc_EUR_trusted
  • effets fixes : year
Summary
# Transformation de la variable année
bso_pop <- bso_pop %>% 
    mutate(year = year - 2013)

# Modèle linéaire
lm1 <- lm(amount_apc_EUR_trusted ~ year, data = bso_pop)

# Summary du modèle
summary <- tidy(lm1)
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
term estimate std.error statistic p.value
(Intercept) 1633.06034 5.168348 315.97340 0
year 39.79084 1.063911 37.40053 0 

# Y estimé 
prediction_lm1 <- round(predict(lm1),2)
    # Mesures d'erreurs
mae_lm1 <- mean(abs(prediction_lm1 - bso_pop$amount_apc_EUR_trusted))
mdae_lm1 <- median(abs(prediction_lm1 - bso_pop$amount_apc_EUR_trusted))
sae_lm1 <- sum(abs(prediction_lm1 - bso_pop$amount_apc_EUR_trusted))

Chaque année, le montant des APC (fiables) augmente de 39.79 euros en moyenne.


Visualisation
require(ggiraph)
require(ggiraphExtra)
require(plyr)
ggPredict(lm1, se=TRUE, interactive=TRUE)


Modèle 2

  • Y : amount_apc_EUR_trusted
  • effets fixes : year + tier + oa_color_finale + is_covered_by_couperin + is_french_CA_bso_wos_openalex_single_lang + bso_classification
# Modèle linéaire
lm2 <- lm(amount_apc_EUR_trusted ~ year + tier + oa_color_finale + is_covered_by_couperin + is_french_CA_bso_wos_openalex_single_lang + bso_classification, data = bso_pop)

# Summary du modèle
summary <- tidy(lm2)
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
term estimate std.error statistic p.value
(Intercept) 1949.00037 11.2763319 172.839925 0.0e+00
year 36.92390 0.9814459 37.621943 0.0e+00
tier2 -271.02121 8.5219874 -31.802583 0.0e+00
tier3 -401.64620 8.4729780 -47.403191 0.0e+00
tier4 -349.58549 9.0370551 -38.683563 0.0e+00
oa_color_finalehybridOA 763.45374 7.1434349 106.874879 0.0e+00
is_covered_by_couperinnon covered -72.01286 8.4917775 -8.480304 0.0e+00
is_french_CA_bso_wos_openalex_single_langnon french CA 71.06395 4.2454420 16.738881 0.0e+00
bso_classificationChemistry -404.66028 12.7847102 -31.651893 0.0e+00
bso_classificationComputer and information sciences -478.30548 16.1292903 -29.654466 0.0e+00
bso_classificationEarth, Ecology, Energy and applied biology -287.59635 7.6822988 -37.436236 0.0e+00
bso_classificationEngineering -587.72553 21.7926464 -26.968984 0.0e+00
bso_classificationHumanities -881.87264 27.4958951 -32.072884 0.0e+00
bso_classificationMathematics -1010.07039 28.5262179 -35.408493 0.0e+00
bso_classificationMedical research -23.60761 4.8664055 -4.851138 1.2e-06
bso_classificationPhysical sciences, Astronomy -404.98640 8.6786828 -46.664500 0.0e+00
bso_classificationSocial sciences -470.37703 27.9288201 -16.841994 0.0e+00 

# Analyse des résidus
lm2_residu <- round(residuals(lm2),2)

# Y estimé 
prediction_lm2 <- round(predict(lm2),2)
    # Mesures d'erreurs
mae_lm2 <- mean(abs(prediction_lm2 - bso_pop$amount_apc_EUR_trusted))
mdae_lm2 <- median(abs(prediction_lm2 - bso_pop$amount_apc_EUR_trusted))
sae_lm2 <- sum(abs(prediction_lm2 - bso_pop$amount_apc_EUR_trusted))


Modèles multi-niveaux : effets fixes + aléatoires


Modèle 1

  • Y : amount_apc_EUR_trusted
  • effets fixes : year
  • effets aléatoires : bso_classification sur intercept

Chaque discipline dispose de son ordonnée à l’origine, et l’année impacte l’angle de la pente de régression.

Summary
## Modèle multi-niveaux
mm1 <- lmer(amount_apc_EUR_trusted ~ year + (1 | bso_classification), data = bso_pop)

# Summary du modèle
print_summary <- function(model, name){
    
    summary <- tidy(model)
    print(kable(summary, format = "html") %>% 
      kable_styling(bootstrap_options = c("striped", "hover"), 
                    full_width = F, 
                    html_font = "sans-serif"))
    # Y estimé 
    prediction <- predict(model)
        # Mesures d'erreurs
    mae <- mean(abs(prediction - bso_pop$amount_apc_EUR_trusted))
    mdae <- median(abs(prediction - bso_pop$amount_apc_EUR_trusted))
    sae <- sum(abs(prediction - bso_pop$amount_apc_EUR_trusted))
        # Sauvegarde dans l'environnement
    assign(glue("summary_{name}"), summary, envir = .GlobalEnv)
    assign(glue("prediction_{name}"), prediction, envir = .GlobalEnv)
    assign(glue("mae_{name}"), mae, envir = .GlobalEnv)
    assign(glue("mdae_{name}"), mdae, envir = .GlobalEnv)
    assign(glue("sae_{name}"), sae, envir = .GlobalEnv)
    
}
print_summary(mm1, "mm1")
effect group term estimate std.error statistic
fixed NA (Intercept) 1305.66655 104.566292 12.48650
fixed NA year 42.21473 1.051654 40.14127
ran_pars bso_classification sd__(Intercept) 329.72755 NA NA
ran_pars Residual sd__Observation 805.72453 NA NA


Visualisation
# Sauvegarde les paramètres du modèle
intercept <- summary_mm1 %>% filter(term == "(Intercept)") %>% select(estimate)
    
# Visualisation des résultats
tidy(mm1, effects = "ran_vals") %>%  # effets aléatoires
    mutate(estimate = intercept$estimate + estimate,  # ajout de l'ordonnée à l'origine
            ymin = estimate - 2 * `std.error`, # intervalles à 95%
            ymax = estimate + 2 * `std.error`) %>% 
    arrange(estimate) %>% 
    mutate(level = as_factor(level)) %>% 
    ggplot(aes(x = level, y = estimate)) +
    ylab("Intercepts estimated") + 
    geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
    geom_hline(yintercept = intercept$estimate, linetype = 2) +
    coord_flip()


Modèle 2

  • Y : amount_apc_EUR_trusted
  • effets fixes : year
  • effets aléatoires : bso_classification sur intercept (ordonnée à l’origine) et trend (pente de la droite de régression)

Chaque discipline dispose de son ordonnée à l’origine, et l’année impacte l’angle de la pente de régression par discipline.

Summary
## Modèle multi-niveaux
mm2 <- lmer(amount_apc_EUR_trusted ~ year + (1 + year | bso_classification), data = bso_pop)

# Summary du modèle
print_summary(mm2, "mm2")
effect group term estimate std.error statistic
fixed NA (Intercept) 1430.33224 101.49955 14.092006
fixed NA year 15.13775 10.39035 1.456905
ran_pars bso_classification sd__(Intercept) 318.00774 NA NA
ran_pars bso_classification cor__(Intercept).year -0.06309 NA NA
ran_pars bso_classification sd__year 31.66700 NA NA
ran_pars Residual sd__Observation 803.31022 NA NA


Visualisation
library(grid)
library(gtable)

viz_intercept_trend <- function(model, name, summary){
    
    # Sauvegarde les paramètres du modèle
    intercept <- summary %>% filter(term == "(Intercept)") %>% select(estimate)
    coef_trend <- summary %>% filter(term == "year") %>% select(estimate)
        # modèle entier
    tidy_model_intercept <- tidy(model, effects = "ran_vals") %>% arrange(estimate) %>% filter(term %in% "(Intercept)") %>% mutate(level = as_factor(level))
    
    # Visualisation des résultats
    
        ## Ordonnées à l'origine
    graph_intercept <- tidy_model_intercept %>% 
        mutate(estimate = intercept$estimate + estimate,  # ajout de l'ordonnée à l'origine
               ymin = estimate - 2 * `std.error`, # intervalles à 95%
               ymax = estimate + 2 * `std.error`) %>% 
        ggplot(aes(x = level, y = estimate)) +
        ylab("Intercepts estimated") + 
        geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
        geom_hline(yintercept = intercept$estimate, linetype = 2) +
        coord_flip()
    
        ## On récupère l'ordre des niveaux pour l'injecter au deuxième graph 
    order_level <- tidy_model_intercept %>% select(level) %>% t()
    
        ## Coefficients des pentes
    graph_trend <- tidy(model, effects = "ran_vals") %>%  # effets aléatoires
        filter(term %in% "year") %>% # pentes
        mutate(level = factor(level, order = TRUE, levels = order_level)) %>% 
        mutate(estimate = coef_trend$estimate + estimate,  # ajout de l'ordonnée à l'origine
               ymin = estimate - 2 * `std.error`, # intervalles à 95%
               ymax = estimate + 2 * `std.error`) %>% 
        ggplot(aes(x = level, y = estimate)) +
        ylab("Slope coefficients estimated") + 
        geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
        geom_hline(yintercept = coef_trend$estimate, linetype = 2) +
        coord_flip() +
        theme(axis.text.y = element_blank(),
              axis.title.y = element_blank())
    
    grid.newpage()
    grid.draw(cbind(ggplotGrob(graph_intercept), ggplotGrob(graph_trend), size="last"))

}
viz_intercept_trend(mm2, "mm2", summary_mm2)


Modèle 3

  • Y : amount_apc_EUR_trusted
  • effets fixes : year
  • effets aléatoires : tier sur intercept (ordonnée à l’origine) et trend (pente de la droite de régression)

Chaque tier dispose de son ordonnée à l’origine, et l’année impacte l’angle de la pente de régression par tier.

Summary
## Modèle multi-niveaux
mm3 <- lmer(amount_apc_EUR_trusted ~ year + (1 + year | tier), data = bso_pop)

# Summary du modèle
print_summary(mm3, "mm3")
effect group term estimate std.error statistic
fixed NA (Intercept) 1821.7423612 324.40976 5.6155597
fixed NA year 23.2517337 25.46493 0.9130884
ran_pars tier sd__(Intercept) 648.7117385 NA NA
ran_pars tier cor__(Intercept).year -0.9085942 NA NA
ran_pars tier sd__year 50.8713498 NA NA
ran_pars Residual sd__Observation 792.1446766 NA NA


Visualisation
viz_intercept_trend(mm3, "mm3", summary_mm3)


Modèle 4

  • Y : amount_apc_EUR_trusted
  • effets fixes : year
  • effets aléatoires : is_covered_by_couperin sur intercept (ordonnée à l’origine) et trend (pente de la droite de régression)

Chaque groupe d’articles selon qu’ils soient couvert par les accords Couperin dispose de son ordonnée à l’origine, et de son coefficient de pente de régression impacté par l’année.

Summary
## Modèle multi-niveaux
mm4 <- lmer(amount_apc_EUR_trusted ~ year + (1 + year | is_covered_by_couperin), data = bso_pop)

# Summary du modèle
print_summary(mm4, "mm4")
effect group term estimate std.error statistic
fixed NA (Intercept) 1887.64180 344.18858 5.484324
fixed NA year 30.58142 10.40395 2.939406
ran_pars is_covered_by_couperin sd__(Intercept) 486.68880 NA NA
ran_pars is_covered_by_couperin cor__(Intercept).year -1.00000 NA NA
ran_pars is_covered_by_couperin sd__year 14.64059 NA NA
ran_pars Residual sd__Observation 792.57366 NA NA


Visualisation
viz_intercept_trend(mm4, "mm4", summary_mm4)


Modèle 5

  • Y : amount_apc_EUR_trusted
  • effets fixes : year
  • effets aléatoires : oa_color_finale sur intercept (ordonnée à l’origine) et trend (pente de la droite de régression)

Chaque couleur OA dispose de son ordonnée à l’origine, et l’année impacte l’angle de la pente de régression par couleur OA.

Summary
## Modèle multi-niveaux
mm5 <- lmer(amount_apc_EUR_trusted ~ year + (1 + year | oa_color_finale), data = bso_pop)

# Summary du modèle
print_summary(mm5, "mm5")
effect group term estimate std.error statistic
fixed NA (Intercept) 1952.36971 438.95248 4.447793
fixed NA year 25.95908 17.97977 1.443794
ran_pars oa_color_finale sd__(Intercept) 620.72870 NA NA
ran_pars oa_color_finale cor__(Intercept).year -1.00000 NA NA
ran_pars oa_color_finale sd__year 25.38926 NA NA
ran_pars Residual sd__Observation 753.27830 NA NA


Visualisation
viz_intercept_trend(mm5, "mm5", summary_mm5)


Modèle 6

  • Y : amount_apc_EUR_trusted
  • effets fixes : year
  • effets aléatoires : is_french_CA_bso_wos_openalex_single_lang sur intercept (ordonnée à l’origine) et trend (pente de la droite de régression)

Chaque groupe d’articles selon que l’auteur correspondant soit français dispose de son ordonnée à l’origine, et de son coefficient de pente de régression impacté par l’année.

Summary
## Modèle multi-niveaux
mm6 <- lmer(amount_apc_EUR_trusted ~ year + (1 + year | is_french_CA_bso_wos_openalex_single_lang), data = bso_pop)

# Summary du modèle
print_summary(mm6, "mm6")
effect group term estimate std.error statistic
fixed NA (Intercept) 1648.06453 164.325970 10.029240
fixed NA year 38.07521 5.932018 6.418592
ran_pars is_french_CA_bso_wos_openalex_single_lang sd__(Intercept) 232.27695 NA NA
ran_pars is_french_CA_bso_wos_openalex_single_lang cor__(Intercept).year -1.00000 NA NA
ran_pars is_french_CA_bso_wos_openalex_single_lang sd__year 8.25383 NA NA
ran_pars Residual sd__Observation 816.34946 NA NA


Visualisation
viz_intercept_trend(mm6, "mm6", summary_mm6)


Comparaison des modèles multi-niveaux


Visualisation des APC estimés par modèle

APC cumulés
# On rassemble les prévisions
fitted <- bso_pop %>% select(year, amount_apc_EUR_trusted) %>% mutate(year = year + 2013)
fitted <- cbind(fitted, prediction_lm1, prediction_lm2, prediction_mm1, prediction_mm2, prediction_mm3, prediction_mm4, prediction_mm5, prediction_mm6)

# Visualisation
graph <- fitted %>% group_by(year) %>% summarise_all(sum) %>% 
    ggplot(aes(x = year)) +
    geom_line(aes(y = amount_apc_EUR_trusted, color = "valeurs observées")) +
    geom_point(aes(y = amount_apc_EUR_trusted), col = "red", size = 1.2) +
    geom_line(aes(y = prediction_lm1, color = "lm1")) +
    geom_line(aes(y = prediction_lm2, color = "lm2")) +
    geom_line(aes(y = prediction_mm1, color = "mm1")) +
    geom_line(aes(y = prediction_mm2, color = "mm2")) +
    geom_line(aes(y = prediction_mm3, color = "mm3")) +
    geom_line(aes(y = prediction_mm4, color = "mm4")) +
    geom_line(aes(y = prediction_mm5, color = "mm5")) +
    geom_line(aes(y = prediction_mm6, color = "mm6")) +
    scale_color_manual(values = c('valeurs observées' = "red",
                                  "lm1" = "#999999",
                                  "lm2" = "#E69F00",
                                  "mm1" = "#56B4E9",
                                  "mm2" = "#009E73",
                                  "mm3" = "#000000",
                                  "mm4" = "#0072B2",
                                  "mm5" = "#D55E00",
                                  "mm6" = "#CC79A7"
                                  )) +
    labs(x = "Année", y = "Montant des APC (fiables)", title = "Somme des APC estimés par les modèles, par année", colour = "Modélisations") +
    scale_y_continuous(labels = scales::comma) +
    scale_x_continuous(breaks=seq(2013, 2020, 1)) +
    theme_classic() +
    theme(legend.position = "right")
ggplotly(graph)


APC moyens
# Visualisation
graph2 <- fitted %>% group_by(year) %>% summarise_all(.funs = c(mean="mean")) %>% mutate(year = as.numeric(as.character(year))) %>% 
    ggplot(aes(x = year)) +
    geom_line(aes(y = amount_apc_EUR_trusted_mean, color = "valeurs observées")) +
    geom_point(aes(y = amount_apc_EUR_trusted_mean), col = "red", size = 1.2) +
    geom_line(aes(y = prediction_lm1_mean, color = "lm1")) +
    geom_line(aes(y = prediction_lm2_mean, color = "lm2")) +
    geom_line(aes(y = prediction_mm1_mean, color = "mm1")) +
    geom_line(aes(y = prediction_mm2_mean, color = "mm2")) +
    geom_line(aes(y = prediction_mm3_mean, color = "mm3")) +
    geom_line(aes(y = prediction_mm4_mean, color = "mm4")) +
    geom_line(aes(y = prediction_mm5_mean, color = "mm5")) +
    geom_line(aes(y = prediction_mm6_mean, color = "mm6")) +
    scale_color_manual(values = c('valeurs observées' = "red",
                                  "lm1" = "#999999",
                                  "lm2" = "#E69F00",
                                  "mm1" = "#56B4E9",
                                  "mm2" = "#009E73",
                                  "mm3" = "#000000",
                                  "mm4" = "#0072B2",
                                  "mm5" = "#D55E00",
                                  "mm6" = "#CC79A7"
                                  )) +
    labs(x = "Année", y = "Montant des APC (fiables)", title = "Moyenne des APC estimés par les modèles, par année", colour = "Modélisations") +
    scale_y_continuous(labels = scales::comma) +
    scale_x_continuous(breaks=seq(2013, 2020, 1)) +
    theme_classic() +
    theme(legend.position = "right")
ggplotly(graph2)

Les graphiques représentant les APC cumulés et moyens montrent tous deux que la meilleure estimation (la courbe s’approchant le plus des valeurs observées en rouge) est le troisième modèle multi-niveaux, “mm3”. Il s’agit du modèle où la source de variation est le tier.


Tableau récaptiulatif

# Tableau récapitulatif
recapitulatif <- data.frame(
    fix.ef = c("year",
               "year + tier + oa_color_finale + is_covered_by_couperin + is_french_CA_bso_wos_openalex_single_lang + bso_classification", 
               "year", 
               "year", 
               "year", 
               "year", 
               "year", 
               "year"),
    ran.ef = c("", 
               "", 
               "bso_classification", 
               "bso_classification", 
               "tier", 
               "is_covered_by_couperin", 
               "oa_color_finale", 
               "is_french_CA_bso_wos_openalex_single_lang"),
    AIC = c(AIC(lm1), AIC(lm2), AIC(mm1), AIC(mm2), AIC(mm3), AIC(mm4), AIC(mm5), AIC(mm6)),
    BIC = c(BIC(lm1), BIC(lm2), BIC(mm1), BIC(mm2), BIC(mm3), BIC(mm4), BIC(mm5), BIC(mm6)),
    MAE = c(mae_lm1, mae_lm2, mae_mm1, mae_mm2, mae_mm3, mae_mm4, mae_mm5, mae_mm6),
    MDAE = c(mdae_lm1, mdae_lm2, mdae_mm1, mdae_mm2, mdae_mm3, mdae_mm4, mdae_mm5, mdae_mm6),
    SAE = c(sae_lm1, sae_lm2, sae_mm1, sae_mm2, sae_mm3, sae_mm4, sae_mm5, sae_mm6))


# Table finale
recapitulatif <- recapitulatif %>% mutate_if(is.numeric, round, digits = 0)
recapitulatif %>% DT::datatable(options = list(searching = F), width = '60%')


Hétérogénéité de la variance

Tests d’hétérogénéité des variances

La variance peut différer d’un groupe à un autre et il convient de le vérifier statistiquement pour le prendre en compte si besoin. Le test de Bartlett ci-dessous teste l’homogénéité de la variance comme hypothèse nulle, nous l’appliquons à la variable de l’année ‘year’ et regardons son résultat :

# Test d'homogénéité de la variance (H0)
output <- bartlett.test(amount_apc_EUR_trusted ~ interaction(year), data = bso_pop)  #hétérogène

# Présentation des résultats
df <- data.frame(value = unlist(output))
tdf <- as.data.frame(t(df))
tdf <- tdf %>% mutate(data.name = str_replace_all(data.name, "_", "\\_"),  #escape special characters
                          data.name = gsub("bso_pop\\$", "", data.name))
formattable(tdf)
statistic.Bartlett’s K-squared parameter.df p.value data.name method
value 1305.81315717348 7 9.19370755033575e-278 amount_apc_EUR_trusted by interaction(year) Bartlett test of homogeneity of variances

La p-value étant très proche de zéro nous rejetons l’hypothèse nulle : la variance est hétérogène. Dans les modèles suivants nous prenons en compte cette hétérogénéité de la variance et les comparons aux modélisations précédentes à l’aide d’une analyse de variance ANOVA.


Modélisations prenant en compte l’hétérogènéité

Nous construisons un modèle en ajoutant des poids (weights) pouvant varier d’une année à l’autre, pour un effet aléatoire par tier d’éditeurs (modèle le plus performant estimé précédemment).

## Modèle multi-niveaux
library(nlme)
hmm1 <- lme(amount_apc_EUR_trusted ~ year, 
            data = bso_pop, 
            random = ~ 1 + year | tier, 
            weights = varIdent(form = ~ 1 | year), 
            control =list(msMaxIter = 1000, msMaxEval = 1000))

# Summary du modèle
print_summary(hmm1, "hmm1")
effect group term estimate std.error df statistic p.value
fixed NA (Intercept) 1813.747058 234.56788 119611 7.732291 0.0000000
fixed NA year 24.371379 19.38816 119611 1.257024 0.2087474
ran_pars tier sd_(Intercept) 469.015194 NA NA NA NA
ran_pars tier cor_year.(Intercept) -0.856127 NA NA NA NA
ran_pars tier sd_year 38.710926 NA NA NA NA
ran_pars Residual sd_Observation 624.928686 NA NA NA NA
var_model varIdent(1 | year) 0 1.000000 NA NA NA NA
var_model varIdent(1 | year) 1 1.080814 NA NA NA NA
var_model varIdent(1 | year) 2 1.292266 NA NA NA NA
var_model varIdent(1 | year) 3 1.263160 NA NA NA NA
var_model varIdent(1 | year) 4 1.333781 NA NA NA NA
var_model varIdent(1 | year) 5 1.376769 NA NA NA NA
var_model varIdent(1 | year) 6 1.324006 NA NA NA NA
var_model varIdent(1 | year) 7 1.240774 NA NA NA NA


Comparaison des modèles : ANOVA

# ANOVA mm3 et hmm1
anova <- anova(mm3, hmm1)

# Présentation des résultats
df <- data.frame(value = unlist(anova))
tdf <- as.data.frame(t(df))
formattable(tdf)
npar Sum Sq Mean Sq F value
value 1 523160.2 523160.2 0.8337304 

#L'analyse de variance sous hypothèse nulle d'égalité de moyennes, montre qu'il n'y a pas de différence significative entre les deux modèles (p-value = ``r round(tdf$`Pr(>F)`, 2)``). Nous conservons donc le modèle initial présenté en partie précédente.


II - Modélisations for french CA articles


# Préparation du jeu de données prêt à modéliser
bso_frenchCA <- bso_mvp %>% 
    filter(!is.na(bso_classification),
           !is.na(oa_color_finale),
           !is.na(is_french_CA_bso_wos_openalex_single_lang),
           !is.na(amount_apc_EUR),
           !is.na(is_covered_by_couperin)) %>% 
    filter(is_french_CA_bso_wos_openalex_single_lang == "french CA", apc_has_been_paid == 1, amount_apc_EUR > 0)

Par rapport à la partie précédente, ici ne sont gardés que les articles ayant un auteur correspondant français. Il y a ainsi 81054 articles avec un CA français, pour lesquels on sait qu’un APC a été payé, et dont la discipline, le tier, la couleur et la couverture par Couperin sont connus (seuls 4 articles avec un CA français et des APC payés ont été mis de côté parce que l’une des variables citées étaient inconnues).


Analyse exploratoire


amount_apc_EUR


Discipline

detach(package:plyr)
table2_y1_class <- bso_frenchCA %>% group_by(year, bso_classification) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR)) %>% ungroup()

  # Plot
graph_y1_class <- ggplot(table2_y1_class, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y1_class)


Tier

table2_y1_tier <- bso_frenchCA %>% group_by(year, tier) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR)) %>% ungroup() %>% mutate(tier = as.character(tier))

  # Plot
graph_y1_tier <- ggplot(table2_y1_tier, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y1_tier)


Couleur OA

table2_y1_color <- bso_frenchCA %>% group_by(year, oa_color_finale) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR)) %>% ungroup() %>% mutate(oa_color_finale = as.character(oa_color_finale))

  # Plot
graph_y1_color <- ggplot(table2_y1_color, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y1_color)


Covered by Couperin

table2_y1_coup <- bso_frenchCA %>% group_by(year, is_covered_by_couperin) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))

  # Plot
graph_y1_coup <- ggplot(table2_y1_coup, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y1_coup)


amount_apc_EUR_trusted


Discipline

table2_y2_class <- bso_frenchCA %>%  group_by(year, bso_classification) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()

  # Plot
graph_y2_class <- ggplot(table2_y2_class, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y2_class)


Tier

table2_y2_tier <- bso_frenchCA %>% group_by(year, tier) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(tier = as.character(tier))

  # Plot
graph_y2_tier <- ggplot(table2_y2_tier, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y2_tier)


Couleur OA

table2_y2_color <- bso_frenchCA %>% group_by(year, oa_color_finale) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(oa_color_finale = as.character(oa_color_finale))

  # Plot
graph_y2_color <- ggplot(table2_y2_color, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y2_color)


Covered by Couperin

table2_y2_coup <- bso_frenchCA %>% group_by(year, is_covered_by_couperin) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))

  # Plot
graph_y2_coup <- ggplot(table2_y2_coup, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph_y2_coup)


Ici encore en prenant comme population seulement les articles avec le CA français, nous décidons de modéliser le montant des APC fiable.

n_old <- nrow(bso_frenchCA)
bso_frenchCA <- bso_frenchCA %>% select(-amount_apc_EUR) %>% filter(!is.na(amount_apc_EUR_trusted))

Le nombre d’articles pour les modélisations est ainsi de 65785. En ne gardant que les sources d’APC fiables, 15269 articles ont été écartés, soit 18.84% des articles avec des APC payés, un CA français et des valeurs connues pour les autres variables.


Dépendances

# Test de corrélation (Spearman car distribution non normale)
output <- cor.test(bso_frenchCA$amount_apc_EUR_trusted, bso_frenchCA$year, method = "spearman") # corrélation significative
df <- data.frame(value = unlist(output))
tdf <- as.data.frame(t(df))
tdf <- tdf %>% mutate(data.name = str_replace_all(data.name, "_", "\\_"),  #escape special characters
                          data.name = gsub("bso_frenchCA\\$", "", data.name))
formattable(tdf)
statistic.S p.value estimate.rho null.value.rho alternative method data.name
value 39956541196575.7 0 0.15791001897175 0 two.sided Spearman’s rank correlation rho amount_apc_EUR_trusted and year 

# Tests de dépendance
    # Année
save_results(chisq.test(bso_frenchCA$year, bso_frenchCA$bso_classification), "1")
save_results(chisq.test(bso_frenchCA$year, bso_frenchCA$tier), "2")
save_results(chisq.test(bso_frenchCA$year, bso_frenchCA$oa_color_finale), "3")
save_results(chisq.test(bso_frenchCA$year, bso_frenchCA$is_covered_by_couperin), "4")
    # Discipline
save_results(chisq.test(bso_frenchCA$bso_classification, bso_frenchCA$tier), "5")
save_results(chisq.test(bso_frenchCA$bso_classification, bso_frenchCA$oa_color_finale), "6")
save_results(chisq.test(bso_frenchCA$bso_classification, bso_frenchCA$is_covered_by_couperin), "7")
    # Tier
save_results(chisq.test(bso_frenchCA$tier, bso_frenchCA$oa_color_finale), "8")
save_results(chisq.test(bso_frenchCA$tier, bso_frenchCA$is_covered_by_couperin), "9")
    # Couleur OA
save_results(chisq.test(bso_frenchCA$oa_color_finale, bso_frenchCA$is_covered_by_couperin), "10")
# Présentation des résultats dans une table
output <- rbind(tdf_1,tdf_2,tdf_3,tdf_4,tdf_5,tdf_6,tdf_7,tdf_8,tdf_9,tdf_10) %>% 
    mutate(p.value = as.numeric(p.value),
           Dependance = case_when(p.value <= 0.05 ~ "Oui",
                                  TRUE ~ "Non"))
# Escape special characters
output <- output %>% mutate(data.name = str_replace_all(data.name, "_", "\\_"),  #escape special characters
                          data.name = gsub("bso_frenchCA\\$", "", data.name))  #remove base name
rownames(output) <- NULL  #remove rownames
# Display table
kable(output, format = "html", caption = "Résultats des tests de dépendance entre les variables qualitatives") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
Résultats des tests de dépendance entre les variables qualitatives
statistic.X-squared parameter.df p.value method data.name Dependance
1179.36582494436 63 0 Pearson’s Chi-squared test year and bso_classification Oui
6674.06953803312 21 0 Pearson’s Chi-squared test year and tier Oui
49.3016149041887 7 0 Pearson’s Chi-squared test year and oa_color_finale Oui
74.7630393039464 7 0 Pearson’s Chi-squared test year and is_covered_by_couperin Oui
4064.54084270166 27 0 Pearson’s Chi-squared test bso_classification and tier Oui
2327.25470251979 9 0 Pearson’s Chi-squared test bso_classification and oa_color_finale Oui
7596.56394938853 9 0 Pearson’s Chi-squared test bso_classification and is_covered_by_couperin Oui
5757.9702867729 3 0 Pearson’s Chi-squared test tier and oa_color_finale Oui
7564.31853293548 3 0 Pearson’s Chi-squared test tier and is_covered_by_couperin Oui
21632.3827056521 1 0 Pearson’s Chi-squared test with Yates’ continuity correction oa_color_finale and is_covered_by_couperin Oui

Ici encore il existe une forte dépendance entre les variables explicatives du modèle, à ne pas inclure dans une même modélisation donc.


Distribution

# Normalité de Y
mean <- mean(bso_frenchCA$amount_apc_EUR_trusted, na.rm = TRUE)   
sd <- sd(bso_frenchCA$amount_apc_EUR_trusted, na.rm = TRUE)  
ggplot(bso_frenchCA, aes(bso_frenchCA$amount_apc_EUR_trusted)) +
      geom_histogram(aes(y=..density..),  color="black", fill = "steelblue", alpha = 0.2) +
    geom_density( color="red", size = 1, adjust = 5) +
    stat_function(fun = dnorm, colour = "Black", size = 1, args = list(mean = mean, sd = sd))  +
    xlim(c(min(bso_frenchCA$amount_apc_EUR_trusted)-1, max(bso_frenchCA$amount_apc_EUR_trusted)+1))+ ggtitle("Distribution du montant des APC fiables") +
  xlab("Montant des APC fiables") + ylab("Densité") +
     theme_classic()

# Test statistique de normalité
out <- JarqueBeraTest(bso_frenchCA$amount_apc_EUR_trusted, robust = FALSE, method = "chisq")  # Y non normal (on rejette H0)
library(formattable)
df <- data.frame(value = unlist(out))
tdf <- as.data.frame(t(df))
formattable(tdf)
statistic.X-squared parameter.df p.value method data.name
value 36052.1176412677 2 0 Jarque Bera Test bso_frenchCA$amount_apc_EUR_trusted

Le montant des APC fiables n’est pas normalement distribué pour les articles ayant un CA français et où des APC ont été payés.

Modèles multi-niveaux

Modèle 1

  • Y : amount_apc_EUR_trusted
  • effets fixes : year
  • effets aléatoires : bso_classification sur intercept (ordonnée à l’origine) et trend (pente de la droite de régression)

Chaque discipline dispose de son ordonnée à l’origine, et l’année impacte l’angle de la pente de régression par discipline.

Summary

## Modèle multi-niveaux
bso_frenchCA <- bso_frenchCA %>% mutate(year = year - 2013)
mm2.1 <- lmer(amount_apc_EUR_trusted ~ year + (1 + year | bso_classification), data = bso_frenchCA)

# Summary du modèle
print_summary(mm2.1, "mm2.1")
effect group term estimate std.error statistic
fixed NA (Intercept) 1372.2979161 90.66522 15.135880
fixed NA year 14.0251387 11.13582 1.259461
ran_pars bso_classification sd__(Intercept) 281.5986411 NA NA
ran_pars bso_classification cor__(Intercept).year 0.0652288 NA NA
ran_pars bso_classification sd__year 33.4758250 NA NA
ran_pars Residual sd__Observation 764.0597813 NA NA


Visualisation

viz_intercept_trend(mm2.1, "mm2.1", summary_mm2.1)


Modèle 2

  • Y : amount_apc_EUR_trusted
  • effets fixes : year
  • effets aléatoires : tier sur intercept (ordonnée à l’origine) et trend (pente de la droite de régression)

Chaque tier dispose de son ordonnée à l’origine, et l’année impacte l’angle de la pente de régression par tier.

Summary

## Modèle multi-niveaux
mm2.2 <- lmer(amount_apc_EUR_trusted ~ year + (1 + year | tier), data = bso_frenchCA)

# Summary du modèle
print_summary(mm2.2, "mm2.2")
effect group term estimate std.error statistic
fixed NA (Intercept) 1732.9556270 212.67692 8.148301
fixed NA year 26.8037745 20.59414 1.301524
ran_pars tier sd__(Intercept) 425.0766712 NA NA
ran_pars tier cor__(Intercept).year -0.9076197 NA NA
ran_pars tier sd__year 41.0629758 NA NA
ran_pars Residual sd__Observation 763.5601474 NA NA


Visualisation

viz_intercept_trend(mm2.2, "mm2.2", summary_mm2.2)


Modèle 3

  • Y : amount_apc_EUR_trusted
  • effets fixes : year
  • effets aléatoires : is_covered_by_couperin sur intercept (ordonnée à l’origine) et trend (pente de la droite de régression)

Chaque groupe d’articles selon qu’ils soient couvert par les accords Couperin dispose de son ordonnée à l’origine, et de son coefficient de pente de régression impacté par l’année.

Summary

## Modèle multi-niveaux
mm2.3 <- lmer(amount_apc_EUR_trusted ~ year + (1 + year | is_covered_by_couperin), data = bso_frenchCA)

# Summary du modèle
print_summary(mm2.3, "mm2.3")
effect group term estimate std.error statistic
fixed NA (Intercept) 1812.33453 360.19889 5.031483
fixed NA year 30.25574 14.24999 2.123210
ran_pars is_covered_by_couperin sd__(Intercept) 509.27854 NA NA
ran_pars is_covered_by_couperin cor__(Intercept).year -1.00000 NA NA
ran_pars is_covered_by_couperin sd__year 20.06310 NA NA
ran_pars Residual sd__Observation 764.77417 NA NA


Visualisation

viz_intercept_trend(mm2.3, "mm2.3", summary_mm2.3)


Modèle 4

  • Y : amount_apc_EUR_trusted
  • effets fixes : year
  • effets aléatoires : oa_color_finale sur intercept (ordonnée à l’origine) et trend (pente de la droite de régression)

Chaque couleur OA dispose de son ordonnée à l’origine, et l’année impacte l’angle de la pente de régression par couleur OA.

Summary

## Modèle multi-niveaux
mm2.4 <- lmer(amount_apc_EUR_trusted ~ year + (1 + year | oa_color_finale), data = bso_frenchCA)

# Summary du modèle
print_summary(mm2.4, "mm2.4")
effect group term estimate std.error statistic
fixed NA (Intercept) 1929.00036 398.59206 4.8395353
fixed NA year 17.50577 24.02843 0.7285438
ran_pars oa_color_finale sd__(Intercept) 563.59894 NA NA
ran_pars oa_color_finale cor__(Intercept).year -1.00000 NA NA
ran_pars oa_color_finale sd__year 33.93122 NA NA
ran_pars Residual sd__Observation 737.07446 NA NA


Visualisation

viz_intercept_trend(mm2.4, "mm2.4", summary_mm2.4)


Modèle 5

  • Y : amount_apc_EUR_trusted
  • effets fixes : year

On fait un simple modèle linéaire.

Summary

## Modèle multi-niveaux
lm2.1 <- lm(amount_apc_EUR_trusted ~ year, data = bso_frenchCA)

# Summary du modèle
print_summary(lm2.1, "lm2.1")
term estimate std.error statistic p.value
(Intercept) 1561.70792 6.506397 240.02654 0
year 41.10933 1.357211 30.28958 0


Comparaison des modèles multi-niveaux

APC cumulés

# On rassemble les prévisions
fitted2 <- bso_frenchCA %>% select(year, amount_apc_EUR_trusted) %>% mutate(year = year + 2013)
fitted2 <- cbind(fitted2, prediction_mm2.1, prediction_mm2.2, prediction_mm2.3, prediction_mm2.4)

# Visualisation
graph2.1 <- fitted2 %>% group_by(year) %>% summarise_all(sum) %>% 
    ggplot(aes(x = year)) +
    geom_line(aes(y = amount_apc_EUR_trusted, color = "valeurs observées")) +
    geom_point(aes(y = amount_apc_EUR_trusted), col = "red", size = 1.2) +
    geom_line(aes(y = prediction_mm2.1, color = "discipline")) +
    geom_line(aes(y = prediction_mm2.2, color = "tier")) +
    geom_line(aes(y = prediction_mm2.3, color = "covered by Couperin")) +
    geom_line(aes(y = prediction_mm2.4, color = "couleur OA")) +
    scale_color_manual(values = c('valeurs observées' = "red",
                                  "discipline" = "#009E73",
                                  "tier" = "#000000",
                                  "covered by Couperin" = "#0072B2",
                                  "couleur OA" = "#D55E00"
                                  )) +
    labs(x = "Année", y = "Montant des APC (fiables)", title = "Somme des APC estimés par les modèles, par année
(Articles avec un auteur correspondant français)", colour = "Modélisations") +
    scale_y_continuous(labels = scales::comma) +
    scale_x_continuous(breaks=seq(2013, 2020, 1)) +
    theme_classic() +
    theme(legend.position = "right")
ggplotly(graph2.1)


APC moyens

# Visualisation
graph2.2 <- fitted2 %>% group_by(year) %>% summarise_all(.funs = c(mean="mean")) %>% mutate(year = as.numeric(as.character(year))) %>% 
    ggplot(aes(x = year)) +
    geom_line(aes(y = amount_apc_EUR_trusted_mean, color = "valeurs observées")) +
    geom_point(aes(y = amount_apc_EUR_trusted_mean), col = "red", size = 1.2) +
    geom_line(aes(y = prediction_mm2.1_mean, color = "discipline")) +
    geom_line(aes(y = prediction_mm2.2_mean, color = "tier")) +
    geom_line(aes(y = prediction_mm2.3_mean, color = "covered by Couperin")) +
    geom_line(aes(y = prediction_mm2.4_mean, color = "couleur OA")) +
    scale_color_manual(values = c('valeurs observées' = "red",
                                  "discipline" = "#009E73",
                                  "tier" = "#000000",
                                  "covered by Couperin" = "#0072B2",
                                  "couleur OA" = "#D55E00"
                                  )) +
    labs(x = "Année", y = "Montant des APC (fiables)", title = "Moyenne des APC estimés par les modèles, par année
(Articles avec un auteur correspondant français)", colour = "Modélisations") +
    scale_y_continuous(labels = scales::comma) +
    scale_x_continuous(breaks=seq(2013, 2020, 1)) +
    theme_classic() +
    theme(legend.position = "right")
ggplotly(graph2.2)

Visuellement parlant, la meilleure estimation semble être celle réalisée à partir du deuxième modèle où la source de variation est le tier.


Erreurs de prévisions par année

error_year <- fitted2 %>% group_by(year) %>% summarise_all(sum) %>% 
    mutate(error_discipline = round(abs(prediction_mm2.1 - amount_apc_EUR_trusted), 0),
           error_tier = round(abs(prediction_mm2.2 - amount_apc_EUR_trusted), 0),
           error_covered_couperin = round(abs(prediction_mm2.3 - amount_apc_EUR_trusted), 0),
           error_OAcolor = round(abs(prediction_mm2.4 - amount_apc_EUR_trusted), 0)) %>% 
    select(1, 7:10) %>% mutate(year = as.character(year)) %>% 
    bind_rows(summarise_all(., ~if(is.numeric(.)) sum(.) else "Somme des erreurs")) %>% 
    bind_rows(summarise_all(., ~if(is.numeric(.)) mean(.) else "Moyenne des erreurs")) %>% 
    mutate_if(is.numeric, round, digits = 0)
error_year %>% DT::datatable(options = list(searching = F), width = '60%')

Ici nous avons calculé le montant total d’APC payés par année selon le BSO (valeurs réelles) et selon les 4 modèles estimés. Nous avons ensuite soustrait les montants réels aux montants prédits par les modèles (en gardant une valeur absolue), de manière à obtenir une erreur de prévision par année. Enfin, nous avons calculé la somme et la moyenne de ces erreurs, ce qui nous permet de constater que le meilleur modèle est celui où le tier est la source de variation / l’effet aléatoire. C’est donc à partir de ce modèle que nous estimons le montant total d’APC pour les années 2021-2030, tendances inchangées.


Tableau récaptiulatif

# Tableau récapitulatif
recapitulatif2 <- data.frame(
    fix.ef = c("year",
               "year", 
               "year", 
               "year"),
    ran.ef = c("bso_classification", 
               "tier", 
               "is_covered_by_couperin", 
               "oa_color_finale"),
   # R2m = c(r.squaredGLMM(mm2.1)[1], r.squaredGLMM(mm2.2)[1], r.squaredGLMM(mm2.3)[1], r.squaredGLMM(mm2.4)[1]),
   # R2c = c(r.squaredGLMM(mm2.1)[2], r.squaredGLMM(mm2.2)[2], r.squaredGLMM(mm2.3)[2], r.squaredGLMM(mm2.4)[2]),
    AIC = c(AIC(mm2.1), AIC(mm2.2), AIC(mm2.3), AIC(mm2.4)),
    BIC = c(BIC(mm2.1), BIC(mm2.2), BIC(mm2.3), BIC(mm2.4)),
    MAE = c(mae_mm2.1, mae_mm2.2, mae_mm2.3, mae_mm2.4),
    MDAE = c(mdae_mm2.1, mdae_mm2.2, mdae_mm2.3, mdae_mm2.4),
    SAE = c(sae_mm2.1, sae_mm2.2, sae_mm2.3, sae_mm2.4))


# Table finale
    # on arrondit les mesures de qualité
recapitulatif2 <- recapitulatif2 %>% mutate_if(is.numeric, round, digits = 0)
    # on affiche le tableau
recapitulatif2 %>% DT::datatable(options = list(searching = F), width = '60%')

On trouve dans cette table récapitulative différentes mesures de la qualité des modèles. Les critères AIC, BIC ainsi que les mesures d’erreurs MAE (mean absolute error), MDAE (median absolute error) et SAE (sum absolute error) doivent être minimisés. Le classement des modèles au vu de ces statistiques n’est pas le même d’une mesure à une autre, ainsi :

  • AIC, BIC, MDAE : selon ces mesures le modèle qui minimise l’erreur est celui avec la couleur OA ;
  • MAE, SAE : selon ces mesures le meilleur modèle est celui avec la discipline.


Tables pivot (drag and drop)


Data des premiers modèles : articles où un APC a été payé

#devtools::install_github(c("ramnathv/htmlwidgets", "smartinsightsfromdata/rpivotTable"))
library(rpivotTable)
bso_pop <- bso_pop %>% mutate(year = year + 2013)
rpivotTable(bso_pop, rows = "year", cols = c("tier"), width = "100%", height = "400px")


Data des deuxièmes modèles : articles où un APC a été payé ET le CA français

bso_frenchCA <- bso_frenchCA %>% mutate(year = year + 2013)
rpivotTable(bso_frenchCA, rows = "year", cols = c("tier"), width = "100%", height = "400px")


III - Prévisions Y : 2021-2030


Max and min increase per year

Data from models (part II)

# Table de min et max par catégorie de variables quali (intégrées aux modèles)

    # Fonction
table_growth <- function(variable, name){
    
    table <- bso_frenchCA %>% group_by(year, {{variable}}) %>% summarise(nb_articles = n()) %>% ungroup() %>% 
        group_by({{variable}}) %>% mutate(diff = nb_articles - lag(nb_articles),
                                          growth = round((nb_articles - lag(nb_articles)) / lag(nb_articles) * 100, 2)) %>% na.omit() %>% 
        arrange(growth) %>% filter(row_number()==1 | row_number()==n()) %>% 
        mutate(mesure = case_when(row_number() == 1 ~ "min",
                                  row_number() == 2 ~ "max")) %>% 
        pivot_wider(names_from = mesure, values_from = c(nb_articles, year, diff, growth)) %>% 
        rename(variable = 1) %>% mutate(variable = as.character(variable))
    
    # Sauvegarde de la table
    assign(glue("growth_{name}"), table, envir = .GlobalEnv)

}

    # Couleur OA
table_growth(oa_color_finale, "color")
    # Tier
table_growth(tier, "tier")
    # Discipline 
table_growth(bso_classification, "discipline")
    # French CA (total car données filtrées uniquement sur is_french_CA == 1)
growth_year <- bso_frenchCA %>% group_by(year) %>% summarise(nb_articles = n()) %>% ungroup() %>% 
    mutate(diff = nb_articles - lag(nb_articles),
           growth = round((nb_articles - lag(nb_articles)) / lag(nb_articles) * 100, 2)) %>% na.omit() %>% 
    arrange(growth) %>% filter(row_number()==1 | row_number()==n()) %>% 
    mutate(mesure = case_when(row_number() == 1 ~ "min",
                              row_number() == 2 ~ "max"),
           variable = "whole data") %>% 
    pivot_wider(names_from = mesure, values_from = c(nb_articles, year, diff, growth)) %>% 
    mutate(variable = as.character(variable))

# Combinaison de ces valeurs dans une seule table
growth_all <- rbind(growth_color, growth_tier, growth_discipline, growth_year) %>% 
    mutate(Variable = case_when(variable == "1" ~ "tier 1",
                                variable == "2" ~ "tier 2",
                                variable == "3" ~ "tier 3",
                                variable == "4" ~ "tier 4",
                                TRUE ~ variable)) %>% 
    rename(`Évolution minimum nombre d'articles` = diff_min,
           `Évolution maximum nombre d'articles` = diff_max,
           `Growth rate minimum` = growth_min,
           `Growth rate maximum` = growth_max,
           `Année du G.R minimum` = year_min,
           `Année du G.R maximum` = year_max) %>% ungroup() %>% 
    select(Variable, `Évolution minimum nombre d'articles`, `Growth rate minimum`, `Année du G.R minimum`, `Évolution maximum nombre d'articles`, `Growth rate maximum`, `Année du G.R maximum`)

growth_all %>% DT::datatable(options = list(pageLength = 17, searching = F))


Whole dataset

# Table de min et max par catégorie de variables quali (intégrées aux modèles)

    # Function
table_growth2 <- function(variable, name){
    
    table <- data %>% group_by(year, {{variable}}) %>% summarise(nb_articles = n()) %>% ungroup() %>% 
        group_by({{variable}}) %>% mutate(diff = nb_articles - lag(nb_articles),
                                          growth = round((nb_articles - lag(nb_articles)) / lag(nb_articles) * 100, 2)) %>% na.omit() %>% 
        arrange(growth) %>% filter(row_number()==1 | row_number()==n()) %>% 
        mutate(mesure = case_when(row_number() == 1 ~ "min",
                                  row_number() == 2 ~ "max")) %>% 
        pivot_wider(names_from = mesure, values_from = c(nb_articles, year, diff, growth)) %>% 
        rename(variable = 1) %>% mutate(variable = as.character(variable))
    
    # Sauvegarde de la table
    assign(glue("growth2_{name}"), table, envir = .GlobalEnv)

}

    # Couleur OA
table_growth2(oa_color_finale, "color")
    # Tier
table_growth2(tier, "tier")
growth2_tier <- growth2_tier %>% arrange(variable)
    # French CA
table_growth2(is_french_CA_bso_wos_openalex_single_lang, "frenchCA")
growth2_frenchCA <- growth2_frenchCA %>% mutate(variable = case_when(variable == "1" ~ "french CA",
                                                                     variable == "0" ~ "non french CA"))
    # Discipline 
table_growth2(bso_classification, "discipline")
    # French CA (total car données filtrées uniquement sur is_french_CA == 1)
growth2_year <- data %>% group_by(year) %>% summarise(nb_articles = n()) %>% ungroup() %>% 
    mutate(diff = nb_articles - lag(nb_articles),
           growth = round((nb_articles - lag(nb_articles)) / lag(nb_articles) * 100, 2)) %>% na.omit() %>% 
    arrange(growth) %>% filter(row_number()==1 | row_number()==n()) %>% 
    mutate(mesure = case_when(row_number() == 1 ~ "min",
                              row_number() == 2 ~ "max"),
           variable = "whole data") %>% 
    pivot_wider(names_from = mesure, values_from = c(nb_articles, year, diff, growth)) %>% 
    mutate(variable = as.character(variable))

# Combinaison de ces valeurs dans une seule table
growth2_all <- rbind(growth2_color, growth2_tier, growth2_frenchCA, growth2_discipline, growth2_year) %>% 
    mutate(Variable = case_when(variable == "1" ~ "tier 1",
                                variable == "2" ~ "tier 2",
                                variable == "3" ~ "tier 3",
                                variable == "4" ~ "tier 4",
                                TRUE ~ variable)) %>% 
    rename(`Évolution minimum nombre d'articles` = diff_min,
           `Évolution maximum nombre d'articles` = diff_max,
           `Growth rate minimum` = growth_min,
           `Growth rate maximum` = growth_max,
           `Année du G.R minimum` = year_min,
           `Année du G.R maximum` = year_max) %>% ungroup() %>% 
    select(Variable, `Évolution minimum nombre d'articles`, `Growth rate minimum`, `Année du G.R minimum`, `Évolution maximum nombre d'articles`, `Growth rate maximum`, `Année du G.R maximum`)

growth2_all %>% DT::datatable(options = list(pageLength = 22, searching = F))


Scénario 2.A : emballement


Un deuxième scénario vise à quantifier le montant des APC en cas “d’emballement” des articles gold : que ce soit en moyenne d’APC ou en nombre d’articles.

Dans un premier temps, cette évolution est quantifiée de la manière suivante :

  • moyenne des APC

    • l’évolution des APC moyens pour les articles gold est multipliée par 2, soit une hausse de 100% (passant de +47.6 à +95.2 euros par an) ;
    • l’évolution des APC moyens pour les articles hybrid reste constante (soit -12.6 euros par an) ;

  • nombre d’articles

    • l’évolution du nombre d’articles gold est multipliée par 1.25, soit une hausse de 25% (passant de +967 à +1209 articles par an) ;
    • l’évolution du nombre d’articles hybrid reste constante (soit 2014 articles, comme en 2020).

Dans un second temps, les prévisions sont calculées en faisant varier l’évolution du nombre d’articles gold ; non plus une hausse de 25%, mais une hausse de 20% puis de 30%.

Pour ces 3 prévisions, la formule pour calculer le montant total des APC est la même que le scénario précédent, à savoir :

montant total = APC moyen x nombre d’articles


Visualisation des prévisions - hausse de 25%

Moyenne APC

plot_color_predictions_2A_1.25 <- function(predictions_gold, predictions_hybridOA, predictions_total,  
                                   predictions_min_gold, predictions_min_hybridOA, predictions_min_total,
                                   predictions_max_gold, predictions_max_hybridOA, predictions_max_total, 
                                   nom_y_label, nom_title){
    
    graph <- predictions_all_color_2A_1.25 %>% ggplot(aes(x=year)) + 
        geom_line(aes(y = {{predictions_gold}}, color = "gold")) +
        geom_line(aes(y = {{predictions_hybridOA}}, color = "hybridOA")) +
        geom_line(aes(y = {{predictions_total}}, color = "whole")) +
        geom_point(aes(y = {{predictions_total}}), col = "red", size = 1.2) +
        geom_ribbon(aes(ymin = {{predictions_min_gold}}, ymax = {{predictions_max_gold}}), alpha = 0.1) +
        geom_ribbon(aes(ymin = {{predictions_min_hybridOA}}, ymax = {{predictions_max_hybridOA}}), alpha = 0.1) +
        geom_ribbon(aes(ymin = {{predictions_min_total}}, ymax = {{predictions_max_total}}), alpha = 0.1) +
        scale_color_manual(values = c("gold" = "#FFCC00",
                                      "hybridOA" = "#0072B2",
                                      "whole" = "red")) +
        theme(legend.position = "right") +
        labs(x = "Année", y = nom_y_label, title = paste("Évolution", nom_title, "par année et par couleur"), color = " ") +
        scale_y_continuous(labels = scales::comma) +
        scale_x_continuous(breaks=seq(2013, 2030, 2)) +
        geom_vline(xintercept = 2020, linetype = "dashed", color = "black", size = .5) +
        theme_classic()
    ggplotly(graph)

}
plot_color_predictions_2A_1.25(mean_APC_gold_2A_1.25, mean_APC_hybridOA_2A_1.25, mean_APC_total, 
                      mean_APC_min_gold_2A_1.25, mean_APC_min_hybridOA_2A_1.25, mean_APC_total_min, 
                      mean_APC_max_gold_2A_1.25, mean_APC_max_hybridOA_2A_1.25, mean_APC_total_max, 
                      "APC moyens", "de la moyenne des APC")


Nombre d’articles

plot_color_predictions_2A_1.25(nb_articles_gold_2A_1.25, nb_articles_hybridOA_2A_1.25, nb_articles_total, 
                      nb_articles_min_gold_2A_1.25, nb_articles_min_hybridOA_2A_1.25, nb_articles_total_min, 
                      nb_articles_max_gold_2A_1.25, nb_articles_max_hybridOA_2A_1.25, nb_articles_total_max, 
                      "Nombre d'articles", "du nombre d'articles")


Montant total

plot_color_predictions_2A_1.25(montant_apc_gold_2A_1.25, montant_apc_hybridOA_2A_1.25, montant_apc_total, 
                      montant_apc_min_gold_2A_1.25, montant_apc_min_hybridOA_2A_1.25, montant_apc_total_min, 
                      montant_apc_max_gold_2A_1.25, montant_apc_max_hybridOA_2A_1.25, montant_apc_total_max, 
                      "Montant total des APC", "du montant total d'APC payés")


Visualisation des prévisions - hausse de 20%

Moyenne APC

plot_color_predictions_2A_1.2 <- function(predictions_gold, predictions_hybridOA, predictions_total,  
                                   predictions_min_gold, predictions_min_hybridOA, predictions_min_total,
                                   predictions_max_gold, predictions_max_hybridOA, predictions_max_total, 
                                   nom_y_label, nom_title){
    
    graph <- predictions_all_color_2A_1.2 %>% ggplot(aes(x=year)) + 
        geom_line(aes(y = {{predictions_gold}}, color = "gold")) +
        geom_line(aes(y = {{predictions_hybridOA}}, color = "hybridOA")) +
        geom_line(aes(y = {{predictions_total}}, color = "whole")) +
        geom_point(aes(y = {{predictions_total}}), col = "red", size = 1.2) +
        geom_ribbon(aes(ymin = {{predictions_min_gold}}, ymax = {{predictions_max_gold}}), alpha = 0.1) +
        geom_ribbon(aes(ymin = {{predictions_min_hybridOA}}, ymax = {{predictions_max_hybridOA}}), alpha = 0.1) +
        geom_ribbon(aes(ymin = {{predictions_min_total}}, ymax = {{predictions_max_total}}), alpha = 0.1) +
        scale_color_manual(values = c("gold" = "#FFCC00",
                                      "hybridOA" = "#0072B2",
                                      "whole" = "red")) +
        theme(legend.position = "right") +
        labs(x = "Année", y = nom_y_label, title = paste("Évolution", nom_title, "par année et par couleur"), color = " ") +
        scale_y_continuous(labels = scales::comma) +
        scale_x_continuous(breaks=seq(2013, 2030, 2)) +
        geom_vline(xintercept = 2020, linetype = "dashed", color = "black", size = .5) +
        theme_classic()
    ggplotly(graph)

}
plot_color_predictions_2A_1.2(mean_APC_gold_2A_1.2, mean_APC_hybridOA_2A_1.2, mean_APC_total, 
                      mean_APC_min_gold_2A_1.2, mean_APC_min_hybridOA_2A_1.2, mean_APC_total_min, 
                      mean_APC_max_gold_2A_1.2, mean_APC_max_hybridOA_2A_1.2, mean_APC_total_max, 
                      "APC moyens", "de la moyenne des APC")


Nombre d’articles

plot_color_predictions_2A_1.2(nb_articles_gold_2A_1.2, nb_articles_hybridOA_2A_1.2, nb_articles_total, 
                      nb_articles_min_gold_2A_1.2, nb_articles_min_hybridOA_2A_1.2, nb_articles_total_min, 
                      nb_articles_max_gold_2A_1.2, nb_articles_max_hybridOA_2A_1.2, nb_articles_total_max, 
                      "Nombre d'articles", "du nombre d'articles")


Montant total

plot_color_predictions_2A_1.2(montant_apc_gold_2A_1.2, montant_apc_hybridOA_2A_1.2, montant_apc_total, 
                      montant_apc_min_gold_2A_1.2, montant_apc_min_hybridOA_2A_1.2, montant_apc_total_min, 
                      montant_apc_max_gold_2A_1.2, montant_apc_max_hybridOA_2A_1.2, montant_apc_total_max, 
                      "Montant total des APC", "du montant total d'APC payés")


Visualisation des prévisions - hausse de 30%

Moyenne APC

plot_color_predictions_2A_1.3 <- function(predictions_gold, predictions_hybridOA, predictions_total,  
                                   predictions_min_gold, predictions_min_hybridOA, predictions_min_total,
                                   predictions_max_gold, predictions_max_hybridOA, predictions_max_total, 
                                   nom_y_label, nom_title){
    
    graph <- predictions_all_color_2A_1.3 %>% ggplot(aes(x=year)) + 
        geom_line(aes(y = {{predictions_gold}}, color = "gold")) +
        geom_line(aes(y = {{predictions_hybridOA}}, color = "hybridOA")) +
        geom_line(aes(y = {{predictions_total}}, color = "whole")) +
        geom_point(aes(y = {{predictions_total}}), col = "red", size = 1.2) +
        geom_ribbon(aes(ymin = {{predictions_min_gold}}, ymax = {{predictions_max_gold}}), alpha = 0.1) +
        geom_ribbon(aes(ymin = {{predictions_min_hybridOA}}, ymax = {{predictions_max_hybridOA}}), alpha = 0.1) +
        geom_ribbon(aes(ymin = {{predictions_min_total}}, ymax = {{predictions_max_total}}), alpha = 0.1) +
        scale_color_manual(values = c("gold" = "#FFCC00",
                                      "hybridOA" = "#0072B2",
                                      "whole" = "red")) +
        theme(legend.position = "right") +
        labs(x = "Année", y = nom_y_label, title = paste("Évolution", nom_title, "par année et par couleur"), color = " ") +
        scale_y_continuous(labels = scales::comma) +
        scale_x_continuous(breaks=seq(2013, 2030, 2)) +
        geom_vline(xintercept = 2020, linetype = "dashed", color = "black", size = .5) +
        theme_classic()
    ggplotly(graph)

}
plot_color_predictions_2A_1.3(mean_APC_gold_2A_1.3, mean_APC_hybridOA_2A_1.3, mean_APC_total, 
                      mean_APC_min_gold_2A_1.3, mean_APC_min_hybridOA_2A_1.3, mean_APC_total_min, 
                      mean_APC_max_gold_2A_1.3, mean_APC_max_hybridOA_2A_1.3, mean_APC_total_max, 
                      "APC moyens", "de la moyenne des APC")


Nombre d’articles

plot_color_predictions_2A_1.3(nb_articles_gold_2A_1.3, nb_articles_hybridOA_2A_1.3, nb_articles_total, 
                      nb_articles_min_gold_2A_1.3, nb_articles_min_hybridOA_2A_1.3, nb_articles_total_min, 
                      nb_articles_max_gold_2A_1.3, nb_articles_max_hybridOA_2A_1.3, nb_articles_total_max, 
                      "Nombre d'articles", "du nombre d'articles")


Montant total

plot_color_predictions_2A_1.3(montant_apc_gold_2A_1.3, montant_apc_hybridOA_2A_1.3, montant_apc_total, 
                      montant_apc_min_gold_2A_1.3, montant_apc_min_hybridOA_2A_1.3, montant_apc_total_min, 
                      montant_apc_max_gold_2A_1.3, montant_apc_max_hybridOA_2A_1.3, montant_apc_total_max, 
                      "Montant total des APC", "du montant total d'APC payés")


Scénario 2.B : apaisement


Un troisième scénario vise à quantifier le montant des APC en cas “d’apaisement”, où la proportion d’articles hybrid diminue tandis que celle des articles gold augmente.

Dans un premier temps, cette évolution est quantifiée de la manière suivante :

  • moyenne des APC

    • l’évolution des APC moyens pour les articles hybrid reste constante (soit -12.6 euros par an) ;
    • l’évolution des APC moyens pour les articles gold reste constante (soit +47.6 euros par an) ;

  • nombre d’articles

    • le nombre d’articles hybrid est multiplié par 0.9 d’une année à l’autre, soit une baisse de 10% (avec 2014 articles hybrid en 2020, il y en aurait alors 1813 en 2021) ;
    • tous les articles qui ne sont pas hybrid sont considérés comme gold.

Dans un second temps, les prévisions sont calculées en faisant varier l’évolution du nombre d’articles gold ; non plus une hausse de 25%, mais une hausse de 20% puis de 30%.


Visualisation des prévisions - baisse de 10%

Moyenne APC

plot_color_predictions_2B_0.9 <- function(predictions_gold, predictions_hybridOA, predictions_total,  
                                   predictions_min_gold, predictions_min_hybridOA, predictions_min_total,
                                   predictions_max_gold, predictions_max_hybridOA, predictions_max_total, 
                                   nom_y_label, nom_title){
    
    graph <- predictions_all_color_2B_0.9 %>% ggplot(aes(x=year)) + 
        geom_line(aes(y = {{predictions_gold}}, color = "gold")) +
        geom_line(aes(y = {{predictions_hybridOA}}, color = "hybridOA")) +
        geom_line(aes(y = {{predictions_total}}, color = "whole")) +
        geom_point(aes(y = {{predictions_total}}), col = "red", size = 1.2) +
        geom_ribbon(aes(ymin = {{predictions_min_gold}}, ymax = {{predictions_max_gold}}), alpha = 0.1) +
        geom_ribbon(aes(ymin = {{predictions_min_hybridOA}}, ymax = {{predictions_max_hybridOA}}), alpha = 0.1) +
        geom_ribbon(aes(ymin = {{predictions_min_total}}, ymax = {{predictions_max_total}}), alpha = 0.1) +
        scale_color_manual(values = c("gold" = "#FFCC00",
                                      "hybridOA" = "#0072B2",
                                      "whole" = "red")) +
        theme(legend.position = "right") +
        labs(x = "Année", y = nom_y_label, title = paste("Évolution", nom_title, "par année et par couleur"), color = " ") +
        scale_y_continuous(labels = scales::comma) +
        scale_x_continuous(breaks=seq(2013, 2030, 2)) +
        geom_vline(xintercept = 2020, linetype = "dashed", color = "black", size = .5) +
        theme_classic()
    ggplotly(graph)

}
plot_color_predictions_2B_0.9(mean_APC_gold_2B_0.9, mean_APC_hybridOA_2B_0.9, mean_APC_total, 
                      mean_APC_min_gold_2B_0.9, mean_APC_min_hybridOA_2B_0.9, mean_APC_total_min, 
                      mean_APC_max_gold_2B_0.9, mean_APC_max_hybridOA_2B_0.9, mean_APC_total_max, 
                      "APC moyens", "de la moyenne des APC")


Nombre d’articles

plot_color_predictions_2B_0.9(nb_articles_gold_2B_0.9, nb_articles_hybridOA_2B_0.9, nb_articles_total, 
                      nb_articles_min_gold_2B_0.9, nb_articles_min_hybridOA_2B_0.9, nb_articles_total_min, 
                      nb_articles_max_gold_2B_0.9, nb_articles_max_hybridOA_2B_0.9, nb_articles_total_max, 
                      "Nombre d'articles", "du nombre d'articles")


Montant total

plot_color_predictions_2B_0.9(montant_apc_gold_2B_0.9, montant_apc_hybridOA_2B_0.9, montant_apc_total, 
                      montant_apc_min_gold_2B_0.9, montant_apc_min_hybridOA_2B_0.9, montant_apc_total_min, 
                      montant_apc_max_gold_2B_0.9, montant_apc_max_hybridOA_2B_0.9, montant_apc_total_max, 
                      "Montant total des APC", "du montant total d'APC payés")


Visualisation des prévisions - baisse de 15%

Moyenne APC

plot_color_predictions_2B_0.85 <- function(predictions_gold, predictions_hybridOA, predictions_total,  
                                   predictions_min_gold, predictions_min_hybridOA, predictions_min_total,
                                   predictions_max_gold, predictions_max_hybridOA, predictions_max_total, 
                                   nom_y_label, nom_title){
    
    graph <- predictions_all_color_2B_0.85 %>% ggplot(aes(x=year)) + 
        geom_line(aes(y = {{predictions_gold}}, color = "gold")) +
        geom_line(aes(y = {{predictions_hybridOA}}, color = "hybridOA")) +
        geom_line(aes(y = {{predictions_total}}, color = "whole")) +
        geom_point(aes(y = {{predictions_total}}), col = "red", size = 1.2) +
        geom_ribbon(aes(ymin = {{predictions_min_gold}}, ymax = {{predictions_max_gold}}), alpha = 0.1) +
        geom_ribbon(aes(ymin = {{predictions_min_hybridOA}}, ymax = {{predictions_max_hybridOA}}), alpha = 0.1) +
        geom_ribbon(aes(ymin = {{predictions_min_total}}, ymax = {{predictions_max_total}}), alpha = 0.1) +
        scale_color_manual(values = c("gold" = "#FFCC00",
                                      "hybridOA" = "#0072B2",
                                      "whole" = "red")) +
        theme(legend.position = "right") +
        labs(x = "Année", y = nom_y_label, title = paste("Évolution", nom_title, "par année et par couleur"), color = " ") +
        scale_y_continuous(labels = scales::comma) +
        scale_x_continuous(breaks=seq(2013, 2030, 2)) +
        geom_vline(xintercept = 2020, linetype = "dashed", color = "black", size = .5) +
        theme_classic()
    ggplotly(graph)

}
plot_color_predictions_2B_0.85(mean_APC_gold_2B_0.85, mean_APC_hybridOA_2B_0.85, mean_APC_total, 
                      mean_APC_min_gold_2B_0.85, mean_APC_min_hybridOA_2B_0.85, mean_APC_total_min, 
                      mean_APC_max_gold_2B_0.85, mean_APC_max_hybridOA_2B_0.85, mean_APC_total_max, 
                      "APC moyens", "de la moyenne des APC")


Nombre d’articles

plot_color_predictions_2B_0.85(nb_articles_gold_2B_0.85, nb_articles_hybridOA_2B_0.85, nb_articles_total, 
                      nb_articles_min_gold_2B_0.85, nb_articles_min_hybridOA_2B_0.85, nb_articles_total_min, 
                      nb_articles_max_gold_2B_0.85, nb_articles_max_hybridOA_2B_0.85, nb_articles_total_max, 
                      "Nombre d'articles", "du nombre d'articles")


Montant total

plot_color_predictions_2B_0.85(montant_apc_gold_2B_0.85, montant_apc_hybridOA_2B_0.85, montant_apc_total, 
                      montant_apc_min_gold_2B_0.85, montant_apc_min_hybridOA_2B_0.85, montant_apc_total_min, 
                      montant_apc_max_gold_2B_0.85, montant_apc_max_hybridOA_2B_0.85, montant_apc_total_max, 
                      "Montant total des APC", "du montant total d'APC payés")


Visualisation des prévisions - baisse de 5%

Moyenne APC

plot_color_predictions_2B_0.95 <- function(predictions_gold, predictions_hybridOA, predictions_total,  
                                   predictions_min_gold, predictions_min_hybridOA, predictions_min_total,
                                   predictions_max_gold, predictions_max_hybridOA, predictions_max_total, 
                                   nom_y_label, nom_title){
    
    graph <- predictions_all_color_2B_0.95 %>% ggplot(aes(x=year)) + 
        geom_line(aes(y = {{predictions_gold}}, color = "gold")) +
        geom_line(aes(y = {{predictions_hybridOA}}, color = "hybridOA")) +
        geom_line(aes(y = {{predictions_total}}, color = "whole")) +
        geom_point(aes(y = {{predictions_total}}), col = "red", size = 1.2) +
        geom_ribbon(aes(ymin = {{predictions_min_gold}}, ymax = {{predictions_max_gold}}), alpha = 0.1) +
        geom_ribbon(aes(ymin = {{predictions_min_hybridOA}}, ymax = {{predictions_max_hybridOA}}), alpha = 0.1) +
        geom_ribbon(aes(ymin = {{predictions_min_total}}, ymax = {{predictions_max_total}}), alpha = 0.1) +
        scale_color_manual(values = c("gold" = "#FFCC00",
                                      "hybridOA" = "#0072B2",
                                      "whole" = "red")) +
        theme(legend.position = "right") +
        labs(x = "Année", y = nom_y_label, title = paste("Évolution", nom_title, "par année et par couleur"), color = " ") +
        scale_y_continuous(labels = scales::comma) +
        scale_x_continuous(breaks=seq(2013, 2030, 2)) +
        geom_vline(xintercept = 2020, linetype = "dashed", color = "black", size = .5) +
        theme_classic()
    ggplotly(graph)

}
plot_color_predictions_2B_0.95(mean_APC_gold_2B_0.95, mean_APC_hybridOA_2B_0.95, mean_APC_total, 
                      mean_APC_min_gold_2B_0.95, mean_APC_min_hybridOA_2B_0.95, mean_APC_total_min, 
                      mean_APC_max_gold_2B_0.95, mean_APC_max_hybridOA_2B_0.95, mean_APC_total_max, 
                      "APC moyens", "de la moyenne des APC")


Nombre d’articles

plot_color_predictions_2B_0.95(nb_articles_gold_2B_0.95, nb_articles_hybridOA_2B_0.95, nb_articles_total, 
                      nb_articles_min_gold_2B_0.95, nb_articles_min_hybridOA_2B_0.95, nb_articles_total_min, 
                      nb_articles_max_gold_2B_0.95, nb_articles_max_hybridOA_2B_0.95, nb_articles_total_max, 
                      "Nombre d'articles", "du nombre d'articles")


Montant total

plot_color_predictions_2B_0.95(montant_apc_gold_2B_0.95, montant_apc_hybridOA_2B_0.95, montant_apc_total, 
                      montant_apc_min_gold_2B_0.95, montant_apc_min_hybridOA_2B_0.95, montant_apc_total_min, 
                      montant_apc_max_gold_2B_0.95, montant_apc_max_hybridOA_2B_0.95, montant_apc_total_max, 
                      "Montant total des APC", "du montant total d'APC payés")


Scénario 3 : full OA


Un dernier scénario vise à quantifier le montant des APC dans une perspective full OA.

Un troisième scénario vise à quantifier le montant des APC dans un scénario “d’apaisement”, où la proportion d’articles hybrid diminue tandis que celle des articles gold augmente.

Cette évolution est quantifiée de la manière suivante :

  • moyenne des APC

    • l’évolution des APC moyens pour les articles hybrid reste constante (soit -12.6 euros par an) ;
    • l’évolution des APC moyens pour les articles gold reste constante (soit +47.6 euros par an) ;

  • nombre d’articles

    • le nombre d’articles hybrid est multiplié par 0.9 d’une année à l’autre, soit une baisse de 10% (avec 2014 articles hybrid en 2020, il y en aurait alors 1813 en 2021) ;
    • tous les articles qui ne sont pas hybrid sont considérés comme gold.


IV - Imputation des valeurs manquantes de Y : 2013-2020


 

Document sous licence ouverte